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I.  SUMMARY 

~  The  Transfer  Function  Model  (TFM)  is  an  extensive  computer  program 
which  is  capable  of  computing  the  downstream  parameters  of  a  jet  engine 
exhaust  plume  for  generalized  initial  conditions,  duct  system  geometry,  and 
liquid  and/or  gas  dilution. 

^For  example,  the  TFM  can  be  used  to  predict  the  opacity  of  a  plume  at 
any  distance  from  the  engine  exhaust  plane,  or  the  change  in  opacity  as  the 
engine  load  is  changed. 

"More  specifically,  the  following  parameters  of  the  exhaust  plume  can 
be  calculated  by  the  TTM  at  any  arbitrary  point  downstream  from  the  jet 
exhaust  plane: 

A/1 )  Particulate  size  distribution, 

t A) Gas  temperature, 

0/3)  Gas  density  (and  water  vapor  density), 

D/u)Gas  velocity^ 

E/$)Droplet  size  distribution  (if  condensation  has  occurred  or  if 
water  has  been  injected), 

F/i,  Light  scattered  and  absorbed  by  each  of  the  various  sized  par¬ 
ticulates. 

Cl  .  ? 

8^7/Total  plume  opacity,^  _ 

The  TFM  can  be  used  to  determine  the  effects  of  changing  any  one  of 
several  control  parameters.  Some  of  the  parameters  which  can  be  modified 
are: 

At  the  exhaust  plane: 

A.  Gas  velocity  (subsonic  or  supersonic) 

B.  Particulate  size  distribution 


C.  Mass  flow  rate 


E.  Water  vapor  density 

Within  the  test  cell: 

A.  Injected  air  (or  steam) 

B.  Injected  liquid  water 

C.  Variable  test-cell  cross  section 

The  cost  of  running  the  TFM  for  a  given  situation  is  generally  less 
than  $100.00.  Thus,  the  TFM  should  prove  to  be  a  very  cost-effective  tool 
for  predicting  jet  test  cell  exhaust  plume  characteristics. 

II.  INTRODUCTION  AND  BACKGROUND 

The  U.S.  Navy  has  a  need  to  predict  the  opacity  of  jet  engine  test 
cell  emissions  and  the  parameters  which  affect  the  test  cell  exhaust  opacity. 

Because  of  the  difficulty  and  cost  in  measuring  the  particulate  size 
and  opacity  of  jet  engine  emission  plumes,  a  computer  code  has  been  devel¬ 
oped  which  predicts  the  opacity.  The  code  incorporates  the  basic  laws  of 
thermodynamics,  hydrodynamics,  and  aerosol  behavior,  together  with  the  laws 
which  govern  light  scattering  and  absorption,  to  describe  both  the  condi¬ 
tions  inside  the  "smoke"  as  it  moves  away  from  the  engine  and  the  opacity  of 
the  plume  itself. 

Since  the  computer  code  describes  the  exhaust  gas  aerosol  from  the 
exhaust  phase  of  the  jet  engine  until  the  aerosol  is  transferred  into  the 
atmosphere,  the  code  is  called  a  Transfer  Function  Model  (TFM). 

In  order  to  make  the  TFM  as  general  as  possible,  allowance  has  been 
made  for  the  injection  of  steam,  air,  other  gases,  or  liquid  water  at  arbi¬ 
trary  points  downstream  from  the  gas.  This  capability  makes  possible  the 
study  of  the  effect  of  such  injections  on  the  opacity  of  the  plume  as  it 
rises  into  the  atmosphere. 


The  original  computer  code  upon  which  the  TFM  is  based  was  a  code  used 
to  simulate  various  pollution  control  devices.  This  aspect  of  the  TFM  would 
make  it  a  valuable  asset  if,  in  the  future,  the  Navy  wanted  to  evaluate  pro¬ 
posed  test  cell  emissions  control  methods.  The  TFM  could  be  used  to  simu¬ 
late  each  of  the  proposed  methods  to  determine  energy  costs,  particulate 
removal  efficiency,  resultant  plume  opacity,  etc. 

Later  sections  of  this  report  describe  in  detail  the  equations  upon 
which  the  TFM  is  based  and  the  computational  techniques  which  are  used.  In 
short,  the  TFM  follows  a  representative  parcel  of  "smoke"  from  the  time  it 
leaves  the  jet  engine  exhaust  until  it  rises  into  the  atmosphere.  As  the 
TFM  "steps"  the  parcel  along,  it  calculates  the  changes  in  the  intrinsic 
parameters  (pressure,  temperature,  etc.)  of  the  gas  and  the  changes  in  the 
size  distribution  and  number  density  of  the  aerosol  particulates.  If  water 
has  been  injected,  the  TFM  calculates  (at  each  step)  the  amount  of  water 
which  evaporates  and  the  size  distribution  and  number  density  of  water  drop¬ 
lets. 

The  second  part  of  th’e  TFM  then  uses  the  aerosol  particulate  and  water 
droplet  size  distributions  and  number  density  to  calculate  the  opacity  of 
the  plume.  The  second  part  of  the  TFM  is  based  on  the  generalized  scatter¬ 
ing  equations  for  electromagnetic  radiation. 

In  developing  the  TFM,  care  was  taken  to  make  it  as  general  as  pos¬ 
sible  without  increasing  the  computer  time  unnecessarily.  The  result  is  a 
model  which  can  be  easily  modified  or  expanded  as  required. 

The  TFM  should  provide  a  very  valuable,  cost-effective  tool  for  use  in 


determ in-'ng  jet  engine  test  facility  emissions. 


III.  THEORETICAL  BASIS  AND  COMPUTATIONAL  TECHNIQUES 


A.  Exhaust  Gas/Aerosol  Dynamics 

The  TFM  calculates  the  intrinsic  parameters  (temperature,  pres¬ 
sure,  density,  velocity,  etc.)  of  the  exhaust  gas  as  it  moves  away  from  a 
jet  engine.  These  calculations  are  based  on  the  basic  principles  of  fluid 
dynamics  and  thermodynamics  for  compressible  fluid  flow  (see  references 
1-4).  Those  basic  principles  include  the  following  laws  and  assumptions: 

1.  Conservation  of  Mass 

2.  Newton's  Second  Law 

3.  First  Law  of  Thermodynamics 

4.  Second  Law  of  Thermodynamics 

5.  Equation  of  State  for  a  Perfect  Gas 

6.  Compressible  Flow 

The  equations  which  represent  these  basic  principles  are  formu¬ 
lated  in  the  method  outlined  by  Shapiro  (reference  1).  Using  this  method,  a 
matrix  of  influence  coefficients  is  formed,  the  elements  of  which  can  be 
used  to  calculate  the  change  in  each  dependent  variable  due  to  changes  in 
the  independent  variables. 

In  this  method,  the  matrix  elements  of  influence  coefficients  are 
first  calculated  using  the  initial  values  of  the  variables  (i.e.,  the  values 
specified  at  the  jet  engine  exhaust  plane).  The  TFM  then  calculates  the 
changes  in  the  independent  variables  between  that  point  and  a  point  lying  a 
short  distance  downstream.  The  matrix  of  influence  coefficients  are  then 
evaluated  for  the  new  point  and  then  used  to  calculate  the  dependent  vari¬ 
ables  at  the  new  point.  This  process  is  repeated  until  the  exhaust  reaches 
a  pre-spec  if ied  point. 
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At  each  of  the  calculation  points  (dependent  variable  evaluation 
points),  several  other  calculations  are  performed.  References  for  the  basis 
of  each  of  these  calculations  can  be  found  in  the  Appendix,  which  consists 
of  a  FORTRAN  listing  of  the  TFM  with  documentation.  They  fall  under  the 
following  categories: 

1.  Liquid  Injection.  The  TFM  allows  for  the  injection  of  liquid 
at  arbitrary  rates  at  arbitrary  points  downstream.  The  TFM  checks  to  see  if 
one  of  these  ports  has  been  passed  in  going  from  one  calculation  point  to 
the  next.  If  it  has,  the  injected  liquid  is  subjected  to  break-up  into 
smaller  droplets  in  a  separate  routine.  The  drag  which  the  injected  liquid 
exerts  on  the  exhaust  gas  is  calculated  as  is  the  cooling  of  the  gas  due  to 
evaporation  of  the  liquid  and  transfer  of  sensible  heat  from  the  gas  to  the 
liquid.  The  terms  which  arise  from  these  calculations  are  added  to  the  cor¬ 
responding  changes  in  the  independent  variables. 

2.  Gas  Injection.  The  TFM  also  allows  for  the  injection  of  any 
gas  (or  gas  mixture)  at  arbitrary  ports  downstream.  At  each  calculation 
point  the  TFM  checks  to  see  if  a  gas  injection  port  has  been  passed,  the  gas 
mass  injection  rate  is  added  to  the  exhaust  gas  mass  flow  rate,  the  momentum 
of  the  injected  gas  in  the  direction  of  the  exhaust  gas  flow  is  added  in, 
and  the  injected  gas  is  assumed  to  be  completely  mixed  with  the  exhaust  gas 
to  determine  the  resultant  specific  heat  capacity  and  temperature.  The  in¬ 
jected  gas  may  have  arbitrary  temperature  and  relative  humidity. 

3.  Aerosol  Agglomeration.  The  size,  distribution,  and  number 
density  of  the  aerosol  particles  can  change  in  several  ways.  At  present, 
the  TFM  assumes  that  they  may  change  as  follow: 
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(1)  Particulate-particulate  collisions 

(2)  Particulate-liquid  droplet  collisions 

(3)  Liquid  droplet-liquid  droplet  collisions 

(4)  Liquid  droplet  evaporation  (or  condensation) 

The  particulate-particulate  collisions  are  assumed  to  occur 
via  random  motion  of  the  particulate  (diffusional  agglomeration). 

The  particulate-liquid  droplet  collisions  are  due  to  both  dif¬ 
fusional  agglomeration  and  the  fact  that  injected  liquid  droplets  may  have 
an  "ordered"  velocity  which  is  different  from  the  exhaust  gas,  in  which  case 
the  collisions  are  due  to  inertial  impaction. 

The  liquid  droplet-liquid  droplet  collisions  are  calculated 
using  both  diffusional  agglomeration  and  inertial  impaction. 

All  of  these  calculations  are  done  at  each  time  interval  cor¬ 
responding  to  each  downstream  calculation  point. 

4.  Liquid  Droplet  Evaporation  (or  Condensation).  The  TFM,  in  its 
present  form,  assumes  that  there  are  no  chemical  reactions  downstream  of  the 
jet  engine  exhaust.  Thus,  it  assumes  that  a  given  chemical  species  is  con¬ 
served.  As  it  now  stands,  the  TM  monitors  only  one  evaporating  or  condens¬ 
ing  species— water.  At  each  computational  point,  the  TFM  calculates  the 
partial  pressure  of  the  water  vapor  present  in  the  gas.  If  that  pressure  is 
greater  than  the  saturation  vapor  pressure  condensation  on  the  particulates, 
then  pre-existing  water  droplets  can  occur;  if  the  existing  vapor  pressure 
is  less  than  the  saturation  vapor  pressure,  then  evaporation  can  occur.  The 
TFM  calculates  the  amount  of  evaporation  (or  condensation)  which  occurs  at 
each  computational  point  for  each  size  class  of  particulates  and  water 
droplets. 


The  net  water  mass  condensed  between  the  two  computational 
points  is  then  calculated  and  the  corresponding  release  of  latent  and  sen¬ 
sible  heat  is  added  to  the  heat  of  the  gas. 

B.  Opacity  and  Electromagnetic  Scattering  Calculations 

The  TFM  can  calculate  the  scattering  (and  absorption)  by  any  arbi¬ 
trary  size  spherical  particle  for  a  given  wavelength  of  electromagnetic 
radiation  and  a  given  angle  of  scatter.  Thus,  by  using  these  equations  to 
calculate  the  extinction  cross-section  for  each  size-class  particle  (par¬ 
ticulate  or  liquid  droplet)  in  the  jet  exhaust  (at  a  given  point  downstream 
from  the  exhaust  plane),  the  total  extinction  cross-section  can  be  ob¬ 
tained.  Thus,  the  opacity  of  the  plume  can  be  found  by  knowing  the  number 
density  of  each  size-class  particle  and  the  diameter  of  the  plume. 

The  equations  upon  which  the  scattering  and  opacity  calculations 
are  based  are  the  solutions  to  the  general  scattering  theory  (see  references 
5-10). 

These  solutions  are  in  the  form  of  a  series  of  expansions,  some  of 
which  are  straightforward  but  others  require  evolved  numerical  techniques  in 
order  to  assure  rapid  convergence. 

The  TFM  must  have  the  following  parameters  for  each  of  the 
scattering  calculations: 

Diameter  of  the  scattering  particle 
2.  Wavelength  of  the  electromagnetic  radiation 
?.  Angle  defined  by  light  source  :  particle  :  observer 
d.  Imaginary  and  real  components  of  the  index  of  refraction  of 
the  particle  or  the  conductivity  and  relative  permittivity  of 


Further  documentation  of  and  references  for  these  calculations  may 


be  found  in  the  Appendix. 

IV.  CONCLUSIONS 

This  model  provides  an  ability  to  predict  test  facility  particulate 
and  opacity  for  any  combination  of  gas  turbine  engine  and  test  facility/ 

control  device  configurations.  The  cost  is  estimated  to  l/?00th  of  the  cost 
to  measure  the  same  parameters. 

The  TFM  is  a  cost-effective  tool  for  the  determination  of  jet  engine 
test  cell  particulate  loading  and  exhaust  opacity.  It  is  both  flexible  and 
general.  The  flexibility  manifests  itself  in  the  ease  with  which  a  specific 
situation  (engine  size  and  load,  test  cell  geometry,  liquid  injection,  etc.) 
can  be  modeled.  It  is  general  in  the  sense  that  other  effects  (such  as 

chemical  reactions,  hydrocarbon  condensation,  etc.)  can  be  added  without 
modifying  the  theory  and  computational  techniques  upon  which  the  TFM  is 
based. 

V.  RECOMMENDATIONS 

Two  types  of  recommendations  are  pertinent  concerning  the  model. 
First,  the  model  must  be  validated  using  actual  engine  exhaust  plane/top  of 
the  exhaust  stack  data.  Validation  will  provide  air  pollution  control 

districts  with  confidence  of  the  model's  accuracy.  Second,  the  following 
additional  efforts  will  make  the  model  easier  to  use  and  more  powerful  when 
evaluating  control  devices. 

A.  Turbulent  Agglomeration 

As  the  TFM  now  exists,  all  agglomeration  of  the  particulates  is 

calculated  using  diffus’onal  agglomeration  equations  only.  However,  turbu¬ 
lent  agql omerati on  should  be  an  important  process  in  the  modification  of  the 
size  distribution  and  number  density,  especially  for  particulates  greater 
than  about  1  micron. 
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B.  Entrainment  of  Air  Into  the  D'ume 

As  the  TFM  now  exists,  the  amount  of  air  entrained  in  the  plume  is 
an  input  parameter.  However,  more  realistic  method  would  be  to  have  the  TFM 
calculate  the  amount  of  air  entrained  in  the  plume  after  it  leaves  the  test 
cell.  The  entrainment  of  air  into  the  plume  and  the  subsequent  dispersion 
is  required  to  predict  opacity. 

C.  Condensation  of  Hydrocarbons 

As  the  TFM  now  stands,  water  is  the  only  liquid/ vapor  substance 
that  is  allowed  to  undergo  phase  transition.  The  evaporation/condensation 
of  any  other  liquid  can  be  added  to  the  TFM  if  the  saturation  vapor  pressure 
of  the  liquid  is  known. 

D.  Chemical  Reactions 

The  TFM  lends  itself  well  to  the  addition  of  the  effects  of 
chemical  reactions.  If  this  feature  were  added,  the  TFM  could  model  the 
effect  of  afterburning  control  devices  on  the  plume  opacity  or  injection  of 
amnonia  for  N0x  control. 
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PROGRAM  BUILD 

DIMENSION  DAK  20) • NAI ( 20  )  z  SLM  <  20  >  / INSLM ( 20 ) 

£  DAK  I)  is  the  initial  diameter  of  the  aerosol  r  articles  in  the 

C  I  th  size  class/  assumed  to  have  < initially)  the  same  temeera- 

C  tu re  as  the  Jet  exhaust. 

C  NAI(I)  is  the  number  density  of  aerosol  particles  initially  in 

C  the  I  th  size  class. 

C  SLM ( I )  is  the  water-soluble  mass  in  each  of  the  aerosol  part- 

C  icles  in  I  th  size  class. 

c  INSLM C I )  is  the  water-insoluble  solid  mass  in  each  of  the  aerosol 

C  particles  in  the  I  th  size  class. 

1 » DWI ( 20 / 20 ) zNWI (20/20) fTUI( 20 ) zNWIT  <  20 ) 

DW.KI/J)  /  NWI(IzJ)  /and  TUIKJl  are  the  diameter /number  dens— 

C  itaz  and  the  temperature/  respectively  of  the  I  th  size  .in  jected 

C  at  the  J  th  port .NWI ( I / J )  can  be  given  J  port/since  it  is  normal— 

C  ized  later. 

NUIIT(I)  is  an  approximation  to  a  typical  water  spray. 

2 / XUP0RT <  20  >  z WMD0T ( 20 ) »  XAP0RT <  20 )  z  AMD0T  <  20 ) / CONSOL ( 20  > 

XWP0RTCJ)  is  the  distance  to  the  J  th  water  injection  port. 

WMDOT(J)  is  the  mass  injection  water  injection  port. 

XAPORT(J)  is  the  distance  to  the  J  th  air  injection  port. 

AhD0T<J)  is  the  mass  injection  rate  of  air  at  the  J  th  port. 

C0NS0L(J)  is  the  concentration  of  soluble  salts  in  the  injected  water* 
3  z  TIAIR<  20 ) zVXAIR(20) zSPHUMD(20) »MASEQW<20) zSPHTAR(20> 

TIAIR(J)  is  the  temperature  of  the  air  injected  at  the  J  th 
port. 

VXAIR < J >  is  the  axial  component  of  the  velocity  of  the  injected  air 
at  the  J  th  port. 

SPHUMD(J)  is  the  specific  humidity  of  the  sas  being  injected  at 
the  J  th  port. It  is  defined  as  the  ratio  of  the  water  vapor  density  to 
the  total  density  of  the  gas?  thus  a  value  of  1.0  corresponds 
to  steam* 

MASEQU(J)  is  the  mass  eauivilant  weisht(or  ks./kmole) 
of  the  gas  .injected  at  the  J  th  port. 

SPHTAR<J)  is  the  specfic  heat(at  constant  pressure)  injected  at  the 
J  th  port/ < Joules/Kg-deS. Kelvin) « 

Now  Dimension  the  working  variable  arrays. 

4 » RA ( 20 ) / RU ( 20 » 20 ) »S0LUMA(20/20>  » INSGLM( 20/20 ) »H20MAS< 20/20) 

RA<I)  is  the  running  radius  of  the  aerosol  particles. 

RbKI/J)  is  the  running  radius  of  the  injected  water  drops. 

SOLUMA(IzJ)  is  the  running  soluble  mass  in  each  of  the(IzJ)  drops* 
INSQLM(IzJ)  is  the  running  insoluble  mass  in  each  of  the(I«J)  drops* 
H20MAS<I/J)  is  the  mass  of  water  in  each  of  thed/J)  drops* 

5/TU<20/20) / 0XU ( 20 / 20 ) / NW ( 20 / 20 ) /NA ( 20 ) / TA < 20 ) 

TW(IzJ)  is  the  running  temperature  of  the  water  drops  which  were 
injected  at  the  J  th  and  which  were  originally  in  the  I  th 
size  class. 

OXW(IzJ)  is  the  running  axial  component  of  the  velocity  of  the  water 
drops  which  were  injected  at  the  J  th  port  and  whicch  were  originally 
in  the  I  th  size  class. 

NW<IzJ>  is  the  running  number  of  water  drops  which  were  originally  of 
the  I  th  size  class  and  were  injected  at  the  J  th  port* (NIKI mat  be 
less  than  NWKI/J)  due  loses  by  collision  or  coagulation  with  larger 
drops  or  particles.) 
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HA <  I  #  J )  i 5  the  running  number  of  aerosol  particles  which  were  orig¬ 
inally  in  the  I  th  size  class. NA<I)  will.#  in  •  v'ene  rsl  -  be  leas  the;.  MAI 
(I)  because  of  loses-  by  collisions  and  coagulation  with  lsrssr 
particles  or  water  drops. 

TA<I)  is  the  running  temperature  of  the  aerosol  particles  #  which  car. 
be  different  from  the  main  das  flow  because  of  condensation  on  the 
particles . 

6  #  KWI J ( 20 )  #  KAI J  ( 20 )  #  FNORM  ( 20  )  »  DRDIF  ( 20 )  ,  UORKAR  <  20 , 20 )  #TRU<20#20> 

KUIJ(J)  is  a  control  array  to  determine  if  the  J  th  water  spray  hoa 
been  injected?  similarily ?  KAIJ(J)  is  for  the  air  injection  ports. 

FNORM(J)  is  a  normalisation  array  for  the  drop  size  distribution 
injected  at  the  J  th  port. 

DRDIF(I)  is  a  typical  drop  diameter  distribution  to  represent 
the  injected  drops. 

UORKAR  ( I  #  J  >  ond  TRW<I#J)  are  temporary  storage  arrays#  used  when 
injecting  hte  water  sprays. 

7  #  WATM ( 20 )  >  TYARDI <  20 ) #  TYNUHD  <  20 ) #  DAPL0T ( 20 ) 

WATM(I)  is  the  array  of  the  condensed  water  mass  on  each  of  the 
I  th  size  aerosol  particles. 

TYARDI  (X)  .is  a  typical  aerosol  diameter  array  used  to  represen 
the  size  distribution  of  the  aerosol  particles. 

TYNUMD(I)  is  an  un-normalized  size  distribution  which#  with 
DAPL0T  is  the  plotting  array  for  the  diameters#  in  microns. 

8#Aa<6) #Bb(6#6) #Cc( A) »DA(20) #DAIPL0<20> 

Aa»Bb#Cc  are  the  Shapiro  matricies  for  the  influence  coefficients. 
DA(I)  is  the  7  running'  diameter  of  the  aerosol  particles. 

DAIPLO  is  the  plotting  array  for  the  initial  diameters#  microns. 

9  »PLTXXX  < 1000 ) #  T1 YYY  < 1000) » V1YYY ( 1000 ) #  RHYYY ( 1000 ) 

REAL  NAI # INSLM , NU I » MASEQU # INSOLM  #  NU  #  NA  #  NU IT  » 
lHeaw#Loadin#Hixl »Mdotl #Mdot2»Msordl »MsQrd2»Molwtl »Molwt2# 

2NRBAR # INSLGN  # M0UINA >  L # MUH20 # NAMAX 

DATA  Lognor  /100# 112.5#32# 13»3#73#0.24 #0.02 #0 #0»0/ 

DATA  DRDIF  /I . #2. #3. #3. #7. » 10 . # 14. #20. #27 . #35 . #47 . #65 . #80 . » 100 . 

1 # 1 40 . »  200 . 1 300 . r 400 ♦ » 700 . » 1000 . / 

DATA  NWIT  /0,l#0.2»0.4»1.5»l.7#2.5#4,0#6.0»7.0# 10.0# 15.0»20.0 
1# 25.0# 35.0# 33.0# 30.0# 25. 0»15.0#1.00#0. 00/ 

DATA  XUP0RT  /0.30 » 2#00 # 3.00»4 . 00 #5 . 00 #6 . 00# 7 .00# 8 .00 # 9. 00# 10 . 00 
1# 11. 00 #12. 00 #13. 00 #14. 00# 0.0# 16. 00 #50. #50. #60. #79./ 

DATA  XAP0RT  /I .50 #2 .30 #3.50#4 . 50»5 .50 #6. 50#7 .50# 8 .50#9 .50# 10.50 
1# 11. 00 #12.00 #13. 00# 14. 00# 15. 00 #15. 00# 17, 00# 58.00 #0.0# 56.00/ 

DATA  WMD0T  /0. »0.00#0.00»0.00»0.00»0.00#0.00#0.00#0.00»0.00 

1 » 0 . 00  #  0 . 00  #  0 . 00  #0.00  #0.00  #0.00  #0.00  #  0 , 00 » 0 . 00  #0.00/ 

DATA  AMD0T  /00.0»0.00#0.00»00.00#0,0#0.00#0.0#0.00#0,00»000.0 
1#0.0»0.0#0.0»0.0#0.0»0.0»0.0»00.0#0.0#0.0/ 

DATA  CONSOL  /0.0»0.0»0.0>0.00»0.0#0.00#0.0»0.0#0.0»0.0 

I»0.0»0.0»0.00f0.0»0.0#0*0»0#0r0#0#0.0#0.0/ 

DATA  TIAIR  /150 . 0# 20 . #20 . #20 . #20 . » 20 . #20 . » 20 . » 20 . » 20 . 

1 #20.0 #20.0 #20.0 #20.0 #20.0 #20.0 #20.0 #20.0 #20.0 #20.0/ 

DATA  VXAIR  /50 . 0 # 0 . 0 #0 , 0»0 . 0»0 . 0 #0 . 0 #0 , 0 » 0. 0# 0 .0 #0. 0 
1#0,0»0.0#0,0»0.0#0.0>0,0»0.0»0.0»0.0»0.0/ 

DATA  SPHUHD  /0 . 10 » 0 . 0 » 0 .0 » 0 ,0 » 0 . 0 # 0 . 0 » 0 , 0 #0 . 0 » 0 . 0 #0 .0 
l#0.0#0.0»0.0#0.0»0.0r0,0#0.0#0.0#0.0»0.0/ 

DATA  MASEQW  /30 . » 30 ,0 » 30 . 0# 30 . 0 » 30 , 0 » 30 , 0 » 30 . 0 # 30 . 0 #30 .0 # 18 , 0 
1 1 30 . » 30 . #  30 , » 30 ♦ #  30 ♦ » 30 . » 30 . » 30 ♦ #30 . » 30  » / 

DATA  SPHTAR  /1006 .» 1006 ,» 100c .# 1006 .» !«06 .» 1006 .# 1006 . 
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1  *  1006. 1 1900 . 1 1006 . r 1006. » 1006 . ? 1006 . ? 1006 . . 1006 .  • 1006 . -  If >06 . 

2. 1006. r 1006. f 1006./ 

DATA  TWI  /20 . 1 20 • 1 20  .  r 20 . • 20 . » 98  »  >  20 . » 20 . *  20 . ?  20 . 
I»20.f20.»20.»20.»20.f20.f20.f20.»20.»20./ 

DATA  TYARDI  /0 . 05 r 0 . 10 1 0 . 20  *  0 . 30 *  0 . 30  *  0 . 75 / 1 . 0  ? 1 . 3 » 1 . 7 > 2 . 3  t 3 . 1 
1 >  4  *  0  *  5 . 0»  7 . 0  f 10 . 0 1 15 . 0  »  20 . 0  >  30 . Or 50 .0 » 1 00 . 0/ 

DATA  TYNUMD  /500 . 1 1000 . ,  500 . . 200 .  »  100 .  >  50  .  » 10 .  •  7  .  ,  4  .  ,  2  .  * 

11'.  »0. 1  *  0.01  *0. 001  t  0.001  »0.0,0.0»0. 0*0. 0.0.0/ 

Now  the  Universal  constants! 

Cva=76A.4  0  specie  heat  at  constant  volume  of  dry  air. 

Raas=8314.34  9  universal  aas  constant. 

Rma=Raas/29.9S  (?  aas  constant/<mol .  wt.  of  stack  sssj. 

Rmu=Raas/18 .015  0  specific  aas  constant  for  water  vapor. 

PI=3. 14159 

L=2.22583E6  0  latent  heat  of  vaporization  of  water. 

Eps-Rma/Rmw 

Gama- 1.37  9  Cp/Cv  for  dry  aas. 

Cpa-Gama*Cv3  0  Cp  for  dry  air. 

Cpw=418A.  9  speefic  heat  for  liouid  water. 

Difl=2.39E-5  9  water  vapor  diffusion  coef.  (.-i2/sec) 

Rhol=l.E3  9  density  of  liouid  water. 

Cond=2.83E-2  9  Thermal  diffusion  coefficient  in  air. 

The  units  of  Cond  are  Joule/<*eter-sec.-ded» Kelvin) 

Beta=0.03A  9  Beta  is  the  condensation  coefficient. 

Sis»=0.072  0  Surface  tension  for  water. 

Emob=2.2E-4  0  Ion  mobility  in  air  at  1  atm. 

Eta=1.66E-5  0  Kinematic  coef.  of  viscosity  for  air  (in2/sec) 

MUH20=1 . 143E-3  0  Viscosity  of  H20  (Ks/meter-sec. > 

ETAH20=1 .143E-A  0  Kinematic  coef.  of  vis.  for  liouid  H20  (m2/sec.) 

TAYLRC— 2.7  0  Taylor's  constant  for  drop  breakup  in  fast  aas  flows. 

B0LTZ=1.38E-23  ©Boltzmann's  constant  (  Joules/dea) 

C0AGC= ( 2 . /3 . ) *B0LTZ*293 . 1 A/Eta 
C0LC0N=<2./9. >*Rhol/Eta 
CfclBRUP=2.5 

CPH20=1980.  0  Specific  heat  at  constant  pressure  for  H20  vapor. 

How  the  initial  intrinsic  parameters  of  the  aerosol. 

Pdens=2.0E3  0  Particles'  mass  density. 

heow=58  0  Eoiuv.  molecular  wt,  of  NaCl. 

Ey=2.5  0  van't  Hoff  Factor* 

Loadin=5.000  0  Particulate  loadina  (Srains/cu  ft.). 

P0=1.005E5  0  Initial  pressure. 

Sphuml=0*10  0  ■  specific  humidity. 

T0=555.  0  Jet  exhaust  temperatijre  in  dess.  Centisrade. 

V0=2600.  0  Jet  exhaust  velocity  in  ft/sec. 

W0=1597.  0  Mass  flow  rate  throush  Jet  enaine?  lbs/sec. 

Xtot=350. 

A0=4 . 0  0  The  initial  cross-sectional  are3  of  the  flow  duct » so . meter . 

Ratio=2.5  0Ratio  is  the  ratio  of  the  max.  to  min.  cross  sect. 

Now  some  of  the  injected  spray  parameters. 

RH0S0L=1 .05E3  0  Density  of  liouid  to  be  injected  ( ks/cu . meter ) 

Calculate  constants  coefficients  and  convert  alL  parameters  to  MKS 
System  of  units. 

P1=P0 

Ac=2.0*Sid/<Rhol#Rmw) 

Bc=3*l8.015/<4*PI*Rhol> 

C3=L*L*Dif  l/<Cond*Rmw*Rmw) 
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C4 -41cPI/fcD:i  r'l^'Rinw 
C 5=4. *P 1/3. 

C  C 1 1  =  0 . 0  i  a .* £ ;n  .1 0 d e ■*  12 . 4 E  i  9 /  ( P d a n s * 4 *PI) 

Cl 2=3/ ( Rho 1*Cpw ) 

C13=L/  <  4-TPI ) 

C14=4#PI TCond 
C 

C  Initialise  pa ricul ate  P3r3raetar  arrays. 

C 

C 

Ul=V0/3 . 2308 
DO  2  J= 1 r 20 

TUI  (J )  =2^3  . lo+TWI  0J  -  ■  . 

TIAIRC  J)=TIAIR<  J>+*273. 16 
TA ( J ) =T 0+273 .16 
2  CONTINUE 

SUMNUM=0 . 0 

PLDMKS=Losdin*2 . 293E-3  0  Convert  paricle  loading  into  Kd/cu . meter 

DO  9  1=1 f 20 

DAI < I )=TYARDI < I ) *1 . E-6 

RA( I )=DAI ( I ) /2  0  Calculate  initial  particle  radii  in  MRS, 

DA ( I ) =DAI ( I ) 

SLM<I>=0. 

UATM ( I ) =0 • 

SMFACT=C5*RA  < I ) **3*Pdena 
INSLM  < I ) =SMFACT-SLM ( I ) 

SUMNUM=SUMNUM+2HFACT  £TYNUMD( I ) 

DO  10  J=t  *  20 
NWI  ( I » J )  =f!UIT  ( I ) 

NW(I»J)=0.0 

DWI ( I  *  J ) =DRDIF ( I ) *1 . OE-6 
RW ( I »  J )  =  ( D W I ( I » J ) /2 . ) 

DMASSE=C5*RW ( I » J ) **3*RH0S0L 
SOLUMA  < I , J  > =C0NS0L <  J ) *DMASSE 
H20MASI I »  J)=DMASSE-SOLUMA ( I » J) 

INSOLM ( I  *  J ) =0 . 

TU(If J)=TUI(J> 

VXU(I| J)=0. 

WORK AR  < I »  J  >  =0 . 

TRW ( I » J  >  =0 ♦ 

10  CONTINUE 

9  CONTINUE 

FACNRM=SUMNUM/PLDMKS 
DO  7  J=1 »  20 
FNORM  <  J ) =0 ♦ 

KWI J  ( J ) =0  • 

KAI J< J)=0. 

DO  8  1  =  1  * 20 

FNORM  <  J ) =FN0RM <  J ) +NWI ( I >  J ) *C5*RH0S0L*RW ( I » J  >  K*3 
8  CONTINUE 

7  CONTINUE 

DO  6  1=1 r 20 

NAI  <  I ) =TYNUMD ( I ) /FACNRM 
NA  ( I ) =NAI < I ) 

6  CONTINUE 

Condw=0 
Time=0 

Dtma::=l .  00E-4 
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J1--5.0E-3 
De  Ltim-Ol 
X=0 
KKK-1. 

Coun t=0 
Lmn-0 

C  Initialise  the  Plot  parameters  J 
DXPLGT-50 . 0 
:<F’L0T=0.0 
KF'LOT =0 
F  L  =  TC  +  27'i  .  1 6 
vJDCTC=WC/2  *  2 
■:  si»RH*esat ;  ri : 

C  «ixl=Eps*El/<Fi-i£l) 

Sphuml  ®Mixl/ ( 1+Mixl )  G  Specific  humidity. 
hixl=Sphuml/ <  I  .  ■-3phu.fi  1 ) 

El=hi::l#F'l/  ( Eps+Mixl ) 

Cv=< 1+1 ,02*Sphuml >*Cva 
Cp=( l+0«9#Sphuml )*Cpa 
Gaml=CP/Cv 

PrVapf  =  ( Eta/Dif  1 )  '*#  <  1 .  /3  .  ) 

BeSin  stepping. 

30  CONTINUE 
KKK=KKK+1 
Lmn=Lmn+l 

Rm=Rina/<  1- <  I-Eps  >*E1/P1 )  0  Gas  constant  for  moist  air. 

IF  (Lmn.GT . 1 )  GO  TO  11 
' Gasdnl=Pl/(Rm*Tl )  @  Gas  mass  density. 

Mdotl=Gasdnl*Vl#AO  0  Gas  mass  flow  rate. 

Mdotl=UD0T0 

Vl=Mdotl/(Gasdnl*A0) 

Csardl=Gaml*Pl/Gasdnl 

Mso»rdl=Vl*Vl/Csardl 

Snspdl=Csardl** . 5  0  Speed  of  sound  in  gas.' 

Rhowl=El/(Rmw*Tl> 

Rhodal=Gasdnl-Rhowl 

Molwtl=Gasdnl/(Rhodal/29 . 98+Rhowl/18 . 015) 

Areal=A0 
11  CONTINUE 

Calculate  factors  for  diffusion  rate  on  small  droplets. 
Frepth=Dif 1/Beta* < 2. *PI/<Rmw*Tl ) >**(0.5) 

Set  values  for  heat  ventilation  factors. 
PrHeat=<Eta*Cp#Gasdril/Cond)#*<  1  ./3. ) 

IF  (X.LT.XPLOT)  GO  TO  3 
XPLOT =XF’LQT +DXPL0T 
KPL0T=>KPL0T+1 

IF  (KPLOT.GT, 1000)  GO  TO  100 
F’LTXXX  ( KPLQT  )  =X 
T1 YYY  <KPL0T ) =T 1-273. 16 
V1YYY  (KPL0T);=U1 
RHYYY ( KPL0T ) =Gasdnl 
PRINT  202 

PRINT  *,  (PLTQUT<MPIf  KF’LOT)  »«PI»1»5> 

PRINT  *r  <PLTOUT<MPI  .KPLQT)  »MPI~A»B> 

PRINT  *f  <  F’LTOUT  <MP.I » KPLQT ) » MF’I~9*  10) 

IF  (KPLOT.GT. 1 )  GO  TO  333 


A-6 


T».V  ’»  V.v?  V *> ,  *> »>  »> - v  v 


^yrryr^rrrr  •■.•■..■ 


m 

M. 


P9C 

<  . 
«\V 


•:;  '- l  (■:.!,  •sea  i  r  l  .  CHJid  in  .j'-cr.  u.H'i  ?  r  .  w..- :  ?s  -,.:-o  vu;  '  r.d  !•;,-. 

DO  40  J--1.20 

lFiX.LT.XUPORK  J)  :■  GO  TO  oO  I?  .©e  if  J  th  .-art  u  ,s 
I F  (  K W I J <  J  )  .  N  E .  0 )  GO  TO  40  @  Check  to  see  if  injected  os,  vr«v  i. 

KUIJ(J>-J  0  Set  this  water  port  flea  non  zero. 

C  Now  inject  J  th  spray  and  shatter  drop-. 

HO  41  1=1 » 20 

NUI ( I f J)=NUI < I . J ) /FNORM ( J  ) 

C  In  order  to  obtain  the  number  density  of  drops  •  n  the  I  ^  h  si. 

C  range »  injected  at  the  J  th  po r t,  r  multiply  i'WI  <  I  •<  J  )  by  the  water 

C  mass  injection  rate  at  the  J  th  port  *  and  divide  by  tin  total 

C  (aerosol  gas  plus  injected  aas  Plus  evaporated  gas  '■  voLu;(.e  i  l  :*»w 

C  rare  of  the  aas  (  cu.  meters/sec. ; 

Nlv  (if  J'  =NWI  ( 1 ?  J  >  *UMD0T  (  J  )  -  (  rtdoti/Phowl  ) 

C  NU  < I f  J ) =NWI ( I f  J ) KUMDOT <  J ) / ( Mdot 1 /Rhou 1 ) 

C  Check  Taylor's  criteria  for  drop  breakup. 

UREL=U1-UXW  < I f  J ) 

TEST=Rhowl*UREL**2/(2.*Sia/RW< I. J) ) 

IF  ( TEST .LT . TAYLRC )  GO  TO  44 
C  Drop  passes  criteria?  break  it  up. 

U0LIJ=C5*RW< I» J)##3*NU< I , J) 

RBAR=CUBRUP*<2.*RU< I f J)  )**<  . 5 ) * ( MUH20* ( Gig/Rhol ) He*  . 5 
l/<Gasdr.l*YREL**2)  )**< 1 ./3. ) 

NRBAR=VOLI J/ ( C5*RBAR**3 ) 

DO  42  K=1 *  19 

TRCE=(RU(Kf J)+RW(K+lf J) )/2. 

IF  (RBAR.GT.TRCE)  GO  TO  42 

JAY=K 

GO  TO  43 

42  CONTINUE 
JAY=20 

43  CONTINUE 

TRW< JAY» J)=( (C5*U0RKAR< JAYf J)*TRU( JAYf J)**3tU0LIJ>/ 

1 < ( WORKAR  <  JA Y » J ) +NRBAR  >  *C5 ))*#<l«/3.) 

UORKAR ( JAY » J ) =W0RKAR < I » J ) +NRBAR 
GO  TO  41 

44  CONTINUE 

TRW(I» J)=< (UORKAR<I»J)*TRU(If J)**3*NW(I»J)*RUCIf J>**3>/ 

1 ( UORKAR ( I f  J ) +NW < I , J ) ) >  ** < 1 . /3 . ) 

WORKAR ( I r  J ) “UORKAR ( 1 1 J ) +NW <  I  f  J  > 

41  CONTINUE 

DO  45  1  =  1 » 20 

NU ( I f  J ) “UORKAR ( I f  J ) 

RU ( I f  J ) =TRU ( I » J ) 

DMASSE=C5*RH0S0L*RW ( I f  J ) **3 
SOLUMA  < I f  J ) =DMASSE*C0NS0L ( J ) 

H20MAS ( I » J ) “DMASSE-SOLUMA < I f  J  > 

45  CONTINUE 
40  CONTINUE 

C  Set  Dt  for  this  iteration. 

C  Dt=AMAXl ( Dtmax f  Dt ) 

C  DO  200  J“1 f  20 

C  IF  <KUIJ<  J) .E0-0)  GO  TO  200 

C  IF  (UMD0T< J) .EG. 0.0)  GO  TO  200 

C  DO  201  I“1 f  20 

C  IF  (NU(IfJ).EQ.O.O)  GO  TO  201 

C  OPRSR=Esat  ( TU ( l  f  J )  >  K 1 1 .  +Ac/RU <  1 1  : 

C  DMDTR=C4*RU  ( I  f  J )  #  (  E 1,  T 1  -vPRSR/TW  <  I .  J 1 
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w  H T R C N -- C  1  4  SR W  (  [, ?  ,1 >  Sc  <  r  W  I.  T  ..  j  !  - T  L  • 

C  I1  T  R  D  1  ~  A  E1  o  (  (  L  #  .0  M  V  I  R  —  H  i  k  A  C  N  )  '  1  ■  ,  b  K  C p  w  4  P ,  ■  o  .  ••  R  J  <  I  < .  ’  '  'f  T. 

C  fit  i.i  r,  j  t  e  P /  Hi  f  R  [i  1 

C  Dt-AMINl (Dt.ntDt) 


01  CONTINUE 
00  CONTINUE 


c 

0e 

i  t  i  m  =  Dt, 

c 

Now  calculate  the 

ccaau 1st i on 

DO 

49  J 1 : - 1 ,20 

IF 

<  KWI J (Jl) . EQ . 0 )  GO 

10  49 

DO 

50  J2=J1,20 

IF 

<  KW I J  <  J2 ) ♦ EQ . 0  >  GO 

V  f  \  CT  , 

c  The  co33uls+  ion  -•  c 'I  3’..!  1  :;te  1  &«■»  •'  ■  I  ■  -.i  •  -  -•  K  •  '  ’ .  •  ,•  '■  >  ;,r.  I  * 

C  where  nl  is  the  number  Jans  its  or"  the  :  n:  !  I  *  r  ;  .*>  L  !  -  »  I.  L 

C  radius  of  the  smaller  p  article »  n 2  and  r2  are  tne  same  paramo, 

C  the  larger  size  particles  »  M  rl , r2)  is  given  by ; 

C  K(  rl  »  r2)»<2/3)*k*t/Et3lt(  <  rl  +  r2>*#2/<  rl #r2 >  >  ,  where  V  is  tn 

C  Boltzmari  constant)  Eta  is  the  viscsity  of  air  at  temperature  T 

DO  51  11=1,20 

IF(NU(I1,J1). EQ . 0 )  GO  TO  51 

DO  52  12=11,20 

IF ( NW (12, J2) •  EQ , 0 )  GO  TO  52 

DN1DT=-C0AGC*<RU< II, J1 )+RW( 12, J2> ) **2/ ( RW ( I 1 , Ji ) *RW ( 12 , J2 ) ) 
1#NU( II , Jl >*NW< 12, J2) 

IF  <RW(I2,J2).LT.RW<I1,J1>>  THEN 
IL0S=I2 
il  =  ip 
I GANE= I 1 
JG= Jl 
ELSE 
IL0S=I 1 
JL  =  Jl 
IGANE=I2 
JG=J2 
END  IF 

ULMC0A=DN1DT*RW(IL0S, JL)**3*C5*Dt 
NW< ILOS, JL)=NU< IL0S, JL)rDNlDT*Dt 

RW ( IGANE , JG )  =  <  < C5*NW ( IGANE , JG ) *RW  < I GANE , JG ) **3-ULMC0A ) / 

1 < NW < IGANE, JGT*C5)>**<l./3.  ) 

SLMGAN=-DNlDT*Dt*S0LUMA< ILOS, JL) 

S0LUMA ( IGANE , JG ) =S0LUMA  < IGANE , JG ) +SLMGAN/NW < IGANE , JG ) 
INSLGN=-DNlDT*Dt*INS0LM(IL0S, JL) 

INS0LM ( IGANE , JG  >  =  IN30Lf1  < IGANE , JG )  +  INSLGN/NU ( IGANE , JG ) 
H20GAN=-DNlDT*Dt*H20MAS ( ILGS , JL ) 

H20MAS  (  IGANE  ,  JG  )  =H20MAS  <  IGANE  ,  JG  )  +H20GAN/NU  <  IGANE,  JG  > 

52  CONTINUE 
51  CONTINUE 

0  Now  coagulate  the  water  drops  and  the  aerosol  particles, 

DO  53  12=1,20 

I F ( NW  < 1 2  »  J  >  »  EQ . 0 )  GO  TO  53 

DO  54  11=1,20 

IF  (NA( II ) . EQ  •  0 )  GO  TO  54 

DNARDT  =-C0AGC*  <.  RA  <  1 1  )  f  PW  ( 1 2 ,  J  )  )  **2/  (  RA  1  1 1  )  *RU  <  T  2  ,  J  )  )  *NA  <  1 1  ) 

1*NU( 12, J) 

U0LAER=nNARDT»RA  (  I  l )  **3*C5*Dt 
NA  < 1 1 )=NA( II ) +DNARDT*Dt 

RW  (  1 2 ,  J  )  =  (  (  C5KNW  (  1 2  ••  j  )  *!;W  (  1 2 ,  J  )  **3-V0LAER  >  '  <  NW  <.  1 2 «  J  '  TCo  >  ; 

L K*<  l . / 3  .  ) 


A-8 


SQL  IMA  (  12-  J  v=SOLUhrt(  1 2  •>  J ONARDf'XIi  i  “S>LM  .  D  >.  MJ  <  X2  >  J) 

INSOLM  < 12  »  J  >  = INSOLM  C  1 2  >  J  •  -DNARDT* Jt* INSLM  <  I  j.  ..'NJ'  :.2*  j) 

H20MAS  <  1 2  *  J  )  -H20MA3  ( 1 2  .  J  >  •-DNARDT  KDt*WATM  <  1 1  ;  /NW  ■  J!  2*  J) 

54  CONTINUE 
53  CONTINUE 
30  CONTINUE 
4?  CONTINUE 

0  Now  cosjS'jiste  the  aerosol  Particles  with  themselves ♦ 

DO  55  11=1 ,20 

IF  <NA< II ) .EQ.O.O)  GO  TO  53 

DO  56  12=11.20 

IF  \  NA  <  1 2  )  .  !IG  .  0  >  GO  TO  56 

DNARD  T=-CUAGC*  <  RA  <  [1 )  +RA  <  12)  '*.;<2/-  a*  A  < II  )  12:  >  r  N  A  :  II  .  r.NAv  12) 

NA  (  1 1  )  =NA(  II)  fDNAREiT  KDt 

RA<  I2)*(  <NA<  12)  *RA  (12)  *K3-DNARD  T*&t.fcRA  (ID  ft3  )  /  >.  NA  (12)?) 

1**< 1 ,/3. ) 

SLM  < 12 ) =SLH (12) -DNARDT#Dt#SLM (ID /NA  < 12 ) 

INSLM< I2)=INSLM( 12 > ~DNARDT*Dt*INSLM ( I 1 >/NA( 12) 

WATM ( 12 ) =UATM  < 12 ) -DNARDT *Dt*WATM< 12  > /NA ( 12 ) 

56  CONTINUE 

55  CONTINUE 
DNTHLA=0 . 

C  DNTHLA.  is  the  energy  added  to  the  gas  (eer  unit  mass  of  gas) 

C  by  the  injected  air  during  Dt. 

DMDQTA=0 . 

C  DMDOTA  is  the  total  mass  rate  of  dry  air  which  is  injected  at 

C  a  given  port. 

M0WINA=0 . 

C  MOUINA  is  the  weighted  molecular  weight  of  the  injected  air 

C  during  Dt, 

DM0MIA=0. 

C  DMOMIA  is  the  total  momentum  (per  momentum  of  the  gas)  which 

C  is  added  to  the  das  by  the  inJecter. 

VAP I N  J=0 , 0 

0  UAF'INJ  is  the  mass  rate  of  water  vapor  injected  during  Dt. 

DO  60  J=t  *  20 

IF  ( X «LT . XAPOR T ( J ) )  GO  TO  60 
IF  < KAI J ( J ) . ME . 0 )  GO  TO  60 

C 

C  Inject  air  at  J  th  port. 

C 

KAI J ( J ) = J 

D  N  T  H  L  A  =  D  N  T  H  L  A  •  i  S  P  H  TAR(J)*(Tl-TIAI  R  (J)  )  (•  V 1*01/2.  )  *  A  MDCjT(J)  /  M  d  o  1 1 
DMDOT  A=DMDQT  A+  < 1 . 0-SPHUMD ( J ) > SAMDQT ( J ) 

OAF IN  J=VAPIN J+SPHUMD  (  J )  *AMD0T  ( J ) 
h I NA=MGW  I NA+AMDOT  <  J )  XcMASEQW  <  J )  *Gasdn  1  /Hdot,  1 
DM0 M T A = D MOM I A r ( 2 . * 0 X A I R ( J ) / V 1 >  * ( A M D 0 T  <  J ) / h d o  1 1 ) 

CUM  fl'MJE 

I  ■  i oalculata  the  collisions  of  injected  wat.:- ;■  drops  with 

' .  =C  s1,  echer  po  't  a  . 

DO  •.  - 1  * 

if  • !  ".j  i..u  j  .•  gf  r"„ 

if  •' umdo  r  •  j  i )  •  £  i  ,•■).;>  ■  :g 

1  MM 

J2~J3  *20 

. '  .J2  :  ET? .  0  )  GO  TO  7  i. 

■  ;  •  :  :  .  >  >  go  to  - 1 
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n  n 


0  A  l  I A  X  —  0  •  ■' 

NAM  A  ■  A 

DU  99  I  =  L  » 20 

IF  (NA( I  1 .EQ .0.0  ■  JO  TO  99 
DA  (  I )  =2  .  *Rh  (  I  > 

DAIPLO  <  I )  -DA  (  I  )  U  ,E6 
NAMAX=AMAX  1  (  NAMAX  ,  NA  (.  I  )  > 

OAMAX=AMAX 1 ( DA MAX  -  DA i I >  , 

99  CONTINUE 

PL0ZZZ=AL0G1 0  ( NAM  AX  >*10. 

NPl.QZZ=INT  ( PLOZZZ+O  .  5  > 

MMMOD=MOD ( NPLOZZ , 5 ) 

NPLOZZ-NPLOZErf -iihMO  D 
PLOZZZ=NPLOZZ/ 10. 

NAMAX=10.**PL0ZZZ 
DAMAX-DAMAX*1 .Eo 
XXSTP=DAMAX/10 . 

YYSTP=NAhAX/10. 

PRINT  *,  NAMAX  ,  DAMAX  , XXSTP  *  YYSTP 
CALL  BGNPL ( KPLOT ) 

CALL  PAGE  (12., 9.0) 

CALL  NOBRDR 

CALL  AREA2D ( 10 . ,7.5) 

CALL  XNAME< 'PARTICULATE  DIAMETER, (Microns) t' ,  100) 

CALL  YNAME ( ' NUMBER  OF  ® ARTICULATES/CUBIC  METER  IN  SIZE  RANGE $ ' 

1,100) 

C  CALL  HEADIN  ('SIZE  DISTRIBUTION  *',100,1.0,1) 

CALL  MESSAG( 'DISTANCE  FROM  ENGINE  IS  *',100,5.5-7.0) 

CALL  REALN0(X,2, 'ABUT' , 'ABUT' ) 

CALL  MESSAG  < '  METERS* ', 100 ,' ABUT ABUT ' ) 

CALL  GRAF (0.0, XXSTP , DAMAX ,0.0, YYSTP , NAMAX ; 

CALL  YLOG (0.0, XXSTP, 0.0, 7. 5) 

CALL  VBARS ( DA , ' BASE ' , NA , 20 ) 

CALL  CURVE  (DAIPLO, NA, 20, 1 ) 

CALL  ENDPL ( KPLO  T ) 

GO  TO  334 

333  CONTINUE 

CALL  BGNPL (KPLOT) 

CALL  PAGE( 12. ,9. ) 

CALL  NOBRDR 

CALL  AREA2IK  10. ,7.5) 

CALL  XNAME( 'PARTICULATE  DIAMETER, ( Microns >*', 100 > 

CALL  YNAME( 'NUMBER  OF  PARTICULATES/CUBIC  METER  IN  SIZE  RANGE* ' 
1,100) 

CALL  HEADIN  ('SIZE  DISTRIBUTION  OF  PARTICULATES* ’,100,1.5,1) 
CALL  MESSAG  ('DISTANCE  FROM  ENGINE  IS  *',100,5.5,7.0) 

CALL  REALNO  ( X , 2 , ' ABUT ' , ' ABUT ' ) 

CALL  MESSAG  ('  METERS* ', 100 ,' ABUT ',' ABUT ' ) 

CALL  GRAF  ( 0 . 0 , XXSTP , DAMAX , 0 . 0 , YYSTP , NAMAX ) 

CALL  MARKER  (0) 

CALL  CURVE  < DAIPLO ,NAI » 20 *  1  * 

CALL  MARKER  <2> 

CALL  CURVE  (  DAFLOT  .  NA  .  20 , 1 
CALL  ENDPL  (KPLOT) 

334  CONTINUE 

IF  ( KPLOT. GEU  GO  TO  100 
3  CONTINUE 
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73 

72 

71 

70 

C 

C 

C 


C 


E  -  0 


EL  12 =R2  d  ■  V-  «  ,  _  : ,  Ji  >-uxw<I2»J2;  ; 


L.i"  .'LUGN *?<u» 1 1  >>  1  ( II f  Jl )  fURELld/..:, 

IF  ( UGi_ ?\R  . LZ . i . 2 1 4  >  fiO  TO  73 
COLEFF=  (  1+0 . 75  KALQG  \  2  ,  .  CQlPAF.  > . '  <  CGLPAR- 1  .21 
t  i  1-7  =•,  REL 1 2:XP  I  *  ( R W  <  U  » J 1 )  W  v  I  i  >  J  i  .*  v  1 2 » 

.  <  i  i  *  J 1  )  -KNlv  i  2  t  J2  .■  kC  jLErY 

i  F  ( R‘W  ( 1 1  ?  Jl  .  L  f  .  Rw  ( 1 2  ?  J  2  .  )  1  HE.  I 


4 j ;  T  .*  <  -  2  > 
J2 )  £i\U>  C  1.2 


I L  ~  I  L 


JL.  =  J 1 
I  (3-. 1 2 
JG=  J2 
ELSE 
IL  =  I2 
JL  =  J2 
IG  =  I  1 


JG=J1 
END  IF 

DELNWl=DNlDT*Dt 

IF  (DELNWI.GT.NWIILfJL) >  DELNU1=NU(IL? JL) 

NW(  ILf JL)=NW<  ILf  JD-DELNUl 
UGLOOL=  C5*DELNU1*RW<IL? JL)**3 

SOLUMA ( IG f JG )  =SOLUMA <  IG f  JG )  +DELNW1  *S0LUi1A ( IL f  JL ) /NW <  IG f  JG ) 

INSOLMCIGf JG)=INSOLM(IGf JG)+DELNW1*INS0LM<IL> JD/NuH IGf JG) 

H20MAS  ( IG  f  JG )  =H20MAS  (  IGf  JG )  +DELNW1*H20MAS  C  IL  f  JL )  /NW  ( IG  f  JG  > 

RW( IGf JG)=( <C5*NW<IGf JG)*RW< IGf JG)**3+C0LV0L>/<C3*NW< IG» JG) ) ) 
l#*<l./3. ) 

CONTINUE 

CONTINUE 

CONTINUE 

CONTINUE 

Now  calculate  the  collisions  between  the  aerosol  particles  and  the 
injected  liauid  drops  due  to  the  difference  in  their  directed  motion? 
us ins  the  same  method  as  above  for  the  collisions  of  water  drops  with 
themselves. 

DO  74  J2=l f 20 

IF  <NWIJ< J2) .EQ.O)  GO  TO  74 
IF  <WMD0T< J2) .EQ.O)  GO  TO  74 
DO  75  12=1 f 20 

IF  <  NW< 12  ?  J2  > ♦ EG . 0 ;  GO  TO  75 
DO  76  1=1 f20 

IF  (NA (I). EG. 0.0)  GO  TO  76 


VF<ELAW=ABS(V1-0XU(  I2f  J2)  ) 
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CGLPAR=CdL"CON*RA  ( I ;  4RA  <  I )  #VRELAW/RU  .12*  J2  ) 

IF  (COLPAR  ,LE, 1 ,214)  GO  TO  7c 

COLEFF=<  1+0 . 7S#AL0G  <  2  ,  f COLPAR  )/ I  COLPAR- 1 .21*1  >  )  **  I  -2  ) 

DNARDT =VRELAW#PI#  < RA  < I > *RA ( I ) +RW < 12 * J2 ) *RU < 12  * J2 ) > #WA < I >  S 
-fat I2» J2)*C0LEFF 
3ELAER*DNARDT*Dt 

:r  ( bELAER .  G T  . NA ( I ) )  DELAER=NA(  I  > 
d.-t  •;  X  )  -MA  (  I )  -BELAER 
C  2 L A Kf 0 = C 5 Mt D E LAER* R A  v  I ;  ■*'  $ 3 

■1...  LUMA  <I2/J2>  -  SOLUM  A  <I2»J2  >  -fBELAESSSLM  ( I )  /MW  <I2*J2> 

.  .-.S  I'LM  1 12  *  J2  )  «I.MSOLM  <  12 1  J2 )  +DSL/sc*4:fc3L.ft  i  I  )  /NW  <  X2r  J2  / 

/"MS'-  12/  J2)=H20MA3(  12  /  22  >  +  BELAERsWATM  <  I )  /M/  :Z2«J2> 

'  J  '.  £2*  J2  •  ~  •  CZ*?V.:  \  12  7  J2  >  uR'j-:  12  -'2  2  •'  C  V  O'  •  CS'-'M,;  ;  ~ ? . ,  ;  •> 


'■■■'  :  J  .'LL 

-/.i  >c- 


C- TM'jfi 

'*/ o *4  calculate  the  evspa  ration *  the  release  of  latent  heat  *  ar.o  the 
cndu.-ticn  of  heat  f ro.it  the  drees  to  the  has  flow* 

The  saij3tiGft-3  used  for  these  c3lc._l3tic.-i-a  are  as  follows! 


!.)  (VAP.PRE3.at  dr  0  3  surf .) =(s3t .use .  at  T  r )  L'i-hV  r-13M/ '  r.fc.f  3  >  2 

2 .  )  ti m/dt  -  <  4*Pi*r*Di  fl/Rmw )  (Eo/To-Er/Tr > 

.  LIdfft/c!t)a<4/3/#PiH:R}-.ol*C»iw*(r*3lt;  >?.'  dTr/dt)  I-  4#?i*Ccnd*r*(Tr-To> 

w;  is  i  :•  r.  ?  1  to  rhi  surface  ton  a.,  on  ot  the  droplet? 

p*  „*  s  3  i-.om  'i.-.-.  :t.  o r  tier  -  th.o  amount  of  soluble  sails  in 

r.i-;e  -ii  :v'.rii  IS  the  mass  rets  o?  sn%l.n  of  a  droMet  of  radius  r* 

Bifl  is  the  water  vapor  diffusion  coefficient*  Rmw  is  the  aas  constant 
■  r  „!3t&r  vapor*  £o  is  the  vapor  pressure  rar  from  the  drop* To  is  the 
temperature  far-  f a n  the  drcp»£r  and  Tr  ere  the  corresponding  ouanti- 
•;.-2s  ,v;  the  surface  of  the  droplet*'.,  is  the  latent  nest  of  ussorisa- 
of  the  .liquid*  Cpw  is  the  specific  heat  capacity  of  the  : louid 
>ho  •/  and  C  :-r>d  La  the  thermal  di  *  fust.  Lon  coefficient  in  air. 


T:-’rCGN-0.  . 

V-ifCQN  is  the  total  amount  of  best  (per  mass  o*'  .das>  conducted 
.via  thermal  diffusion)  from  ins  aerosol  particles  and  in.iected 
water  dr.-,p*  t->  the  3as  durinS  the  time. 

T.l;/rC  -1  :•  s  total  non  t.u.r-  -.  sp  r  unit  momentum  nf  «»as 
which  is  added  to  th-»  ass  bu  evaporation  of  liouid  during 
2-..  -'Shapiro's  2 .  *ssub  l-f  <  dwsucL.-'w  >  > 

-  A  3=0  , 

TCMMA3  is  the  total  mass  condensed  (per  i.nit  ma«s  nf  the 
L--.S  Or, . 

"riUPEN™  0, 

VE’fPEM  is  too  ener«H  (car  uni*,  moss  of  the  dta«)  »dfied  to  U'*>  '•.->« 
' f.he  evapcra-cad  liquid  r  ’.  -s  ?t . 

20  30  1-1*20 

IF  (NA(I) .EO.C.C;  GO  TO  30 

'PRSR=Esat  <TA(I))#(1*  i-Ac/RA  (  J  ;  -C^.s  3c‘*SL.‘1 ;  l  >  /  :  .t-\  •  "!  «  >  ' 

“''(*C*-'CRA  ( I  '  .X  t  Ei/Tl-VPRSR/TA  ( I  > ) 

.■'-.j-ifv  the  evapo  r  ?t.ion  rate  for  small  droplets! 

-  I-  fRwBMBTR.t  (  RA  (  X  >  /  <  RA  *  X  j  +F  re?»  th  >  > 

;  -hiCi'i  ™l.)  '••I’fRA  (  T  V  'K f  i’A  (  f  >  -"T .!.  ) 


BEST 

AVAILABLE  COPY 


A» 


CVi  i  * '  '  ' 

IF  CEVAMrS 
„  ..  \  I'M 

UihTM  (  .[  :  -  , 

RAC  I  )=DAI  '  T  :■  •  .A 
DMDTR=CONDNM/D  l. 

ELSE 

WATMC I )=WATM< I j+CONDNH 

R  A  <  I  >  =  <  <  SLM  <  I  >  +  INSLM  <  I  >  <  b-,  .  C c  SRbo l  '  •'  "  ■  ’  ' 

END  IF 

TCNMAS=TCNMAS+CONDNM*NA  ( I  )  /G-sdA  L 
THTCON=THTCON+HTRACN*NA<  I  >  *Dt/Gasd.-.l 

TEVPEN»TEVPEN-CQNDNM* ( CPH20* < TA C I ) -T 1 )  -L  >.TNA<  I  -'3a  sen;. 

DTRDT  =  ( L*DMDTR-HTRACN  )  /  ( C5*Rhol*Cpw*KA  (  I )  * *3 ) 

TA ( I ) =  TA ( I ) +BTRDT*Dt 
C  IF  ( TA  < I ) . LT .  Dtemp;-: )  TA<  I  )=Dteit.PX 

80  CONTINUE 

DO  81  J=1 ,20 

IF  ( KWI J <  J ) . EO , 0 )  GO  TO  31 
IF  ( WNDOT  <  J ) . EO . 0 . )  GO  TO  81 
DO  83  1=1,20 

IF  <NU(I,J) .EQ.O.O)  GO  TO  33 

VPRSR=Esat  <  TW  <  I » *J  >  >  * <  1  ♦  +Ac/RW ( I ,  J )  -Ey#Be*SOLUMA  ( I ,  J ) 

1/(RU<I, J>**3> ) 

DMDTR=C4*RW<I, J ) * < El/Tl-VPRSR/TWC I , J) ) 

C  Modify  the  evaporation  rate  for  small  droplets  and  ventilation* 

SaRey=< 2 . *RU ( I , J ) *ABS <  VI -VXU C I , J >  5 /Eta  >**(0.5) 
Frossl=1.0+0.276*SaRey*PrVapf 

DMDTR=DMDTR*(RW< I, J)/(RW<I, J)+Freeth> >*Frossl 
HTRACN=C 1 4*RW <  I ,  J )  *  <  TW  ( I ,  J )  -T1  > 

C  Modify  the  heat  flow  rate  caused  by  ventilation? 

Fross2=l .0+0.276#SaRey#PrHeat 
HTRACN=HTRACN*F  ross2 
CONDNM=DMDTR#Dt 
EVAFMS=-CONDNM 

IF  (EVAPMS«GT .H20MASC I , J) )  THEN 
C0NDNM=— H20MAS  < I , J) 

EVAPMS=-CONDNM 
H20MAS ( I , J ) =0 « 0 

RW ( I ,  J >  =  <  CSOLUMACI,  J > + INSOLM < I , J >  ) / ( Pdens*C5 )  >•**(  1  ./3.  ) 

DMDTR=CONDNM/Dt 

ELSE 

H20MAS ( I , J ) =H20MAS ( I , J ) +CONDNM 

RW(I, J)  =  ( ( <  SQLUMAC I , J  >  +  INSOLM  < I , J ) ) /Pdens+H20MAS ( I , J ) /Rhol )  /C5  > 
1*#< 1 ./3. ) 

END  IF 

TCNMAS=TCNMAS+CONDNM*NW< I , J)/Gasdnl 
THTCON=THTCON+HTRACN*NW< I, J)*Dt/Gasdnl 
TEVPEN=TEVPEN-C0NDNM*(CPH20*<TW< I , J)-T1 )-L+ 

1 <V1*V1-VXU< I, J)*VXW< I, J) )/2. > *NW < I , J > /Gasdn 1 
TEVM0M=TEVM0M+ ( 2 , *VXW ( I , J  > /VI ) * ( EVAPMS*NW C I , J ) /Gasdn 1 ) 
DTRDT=<L*DMDTR-HTRACN)/<C5*CPw*Rhol*RW( I , J) **3; 

TU(I, J)=TU< I , J)+DTRDT*Bt 
C  IF  ( TW <  I ,  J )  . LT  .  Dtempx )  TW(  I ,  J )  =Dtemp>: 

33  CONTINUE 

81  CONTINUE 
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~  rr?  ■>  v 


T  7.  V. 


•  ■ ■  - -1 


7  L' : f '■  M  S = — T  C  N  M  A  S 

C  -riVPM3  is  tiie  total  Uauid  mass  •  ••■it  as-  mass  *  r  ■■  r -■  ■ 

C  durmsa  Dt. 

C 

C  Now  calculate  and  collect  all  of  the  change  >  tn  the  i.nder-„-nd.- n  f 

C  variables  for  input  into  the  Shapiro  matri::. 

C 

C  First  calculate  the  work  done  on  draaains  injected  water, 

SHF'DWX=0 . 0 

C  SHF DUX  is  the  total  work  done  bs  the  sas  (per  unit  mass  of  <2 as  > 

C  on  the  injected  water  via  drad  forces  during  Dt. 

DXAKPM=0 . 0 

C  DXAKPM  is  the  total  draa  force  everted  by  t aas  cn  the 

C  injected  water  during  Dt . 

DO  90  J- i ,  20 

IF  (KUIJ( J) .EQ.O)  GO  TO  90 
IF  (UMDOTC  J)  .ELI,  0.0)  00  TO  90 
DO  91  1=1,20 

IF  <NW< I , J) .EG. 0.0)  GO  TO  91 
F0RSI J=DRAGF ( RU( I , J ) , VI , VXU ( I , J ) »Gasdnl ) 
SHPDWX=SHPDUX+FORSIJ*VXtd<  I ,  J)*Dt*NW<  I ,  Ji/Gasdnl 
DXAKPM=DXAKPM+FORSIJ 

DMASSE=INSOLM ( I , J ) +SQLUMA ( I , J ) +H2CMAS ( I , J  > 

VXW( I , J)=UXW< I , J ) +F0RSI J*Dt/DMASSE 
91  CONTINUE 
90  CONTINUE 

DXAKPH=DXAKPM/  ( A real #Gasdnl *01*01/ 2 . ) 

C  Now  sum  Shapiro's  energy  terms. 

DQDUIDH=THTCON-SHPDWX+  ( DNTHLA+TEOPEN  > 

C  Ne;:  collect  Shapiro's  momentum  terms. 

SHFMOM=DXAKFM- ( DMOMI A+TEOMOM ) 

C  Sum  injected  and  evaporated  mass. 

SHAPDW=  (OAPINJ+DMDOTA)  /Mdotl  +  TEOF’flS 
C  Calculate  the  change  in  the  density  and  molecular  wt. 

DRGSIN=DMDOTA*Gasdnl/Mdotl 
INDENS=OAPINJ*Gasdnl/Mdotl 
EVDENS=TE0PMS*Gasdnl 
Gasdri2=Gasdnl+EVDENS+INDENS+DRGSlN 

Molwt2=<  Gasdnl#Molwtl+EVDENS*18 .013+M0WINA ) /Gasdn2 

Rhow2=Rhowl+E0DENS+INDENS 

Rhoda2=Gasdn2-Rhow2 

Mdot2=Mdotl+SHAPDU*Mdotl 

C  Now  set  up  the  matrix  coefficients  and  calculate  the  independ- 

C  ent  variables.  Calculate  deltas  of  the  independent  variables  to 

C  Set  to  state  2. 

C 

Condw=Condw+T  deltm 
Time  =Time+Deltim 
Delx=01*Deltim 
X=X+Delx 

IF  (X.GT.Xtot)  GO  TO  100 
Area2=FNAREA(X, A0, Ratio; 

W2=Rhow2/Rhoda2 

U2=U2/(1.+U2> 

Gam2=<  1  .+0,90*1)2 )/<  1 .  il . 02*U2 ) *Gams 
C 

C  Calculate  the  independent  variables. 

C 
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Cc  <  2  >  - D >J D U  D  H ,  (  \  v  t  T 1  ) 

Cc  (  3 )  =SHF'MOM 
Cc  (  4 )  =SHAF'DW 

Cc ( 5  > -2 • *  <  Holwt2-Mo  1  ut  1 )  /  ( Mol  wt  1  -riio  lwt2) 

Cc<6>=2.*< Gam2-Gam l ) / ( Gai?.2+Gaii.  I  ) 

A7=Msc*rdl 

B7--Gaml 

C  7 = B  7  :*  A  7 

D7=l ,-A7 

E7=B7-1 . 

F  7  =  1 . +C7 
07*1  ,  +  £•'•!( i\~> 

H  7  - 1  *  -  C  7 

Bb< 1 > 1  >=-2.*67/D7 
Bb< 1 i 2 ) =F7/D7 
Bb< 1 »3>=C7*G7/D7 
Bb ( 1 f  4 ) =2 . *F7*G7/D7 
Bb( 1 t 5>*-F7/D7 
Bb  < 1 r  6 ) =- 1  . 

Bb(2f 1 )=-l ,/D7 
Bb  <  2 1 2 )  =  1 ♦ /D7 
Bb(2»3)=C7/(2.*D7) 

Bb ( 2 » 4 ) =F7/D7 

Bb(2»5)=-1 . /D7 

Bb  <  2  *  A) =0 , 0 

Bb<3, 1  )=E7*A7./<2.*D7> 

Bb ( 3 » 2 ) =H7/ ( 2 . *07 ) 

Bb  <  3 1 3 ) =-B7*E7*A7*A7/ <  4 . *D7  > 

Bb ( 3 » 4 ) =-E7*A7*F7/ < 2 . *D7  > 
Bb(3»5)=-H7/(D7*2. ) 

Bb(3rA)=0«5 

Bb(4»l)=E7/D7*A7 

Bb(4,2)=H7/D7 

Bb(4r3)=-B7*E7*A7/(2.*D7) 

Bb(4,4>=-E7*A7*F7/D7 

Bb(4r5)=E7*A7/D7 

Bb<4»6)=0.0 

Bb<5» 1 )=A7/D7 

Bb(5»2)=-1 ,/D7 

Bb  <  5  >  3 )--C7/ <  2 . *D7 ) 

Bb(5?4)=-(B7+1. )*A7/D7 
Bb(5»5)=l «/D7 
Bb ( 5  >  A ) =0 . 0 
Bb<  A» 1 )=C7/B7 
Bb(6»2)=-C7/D7 

Bb ( 6  r 3 ) =-C7* ( 1 , +E7*A7 ) / ( 2 . *D7 ) 

Bb<6»4)=-2**C7*G7/D7 

Bb ( A » 5 ) =C7/D7 

Bb  <  A  >  A  >  =0 . 0 

DO  20  1  =  1. » 6 

Aa ( T ) -0 . 0 

DO  21  J= 1 » A 

Aa  ( I )  =  Bb  (  I  ■>  J  )  *Cc  :  J)  +As  (  I ) 

CONTINUE 

CONTINUE 


Calculate  the  new  valuer;  of  the  intrinsic  parameters  o 
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M«c»rd2=M-3«i'dl  X  (  1  .  •  1  >  ' 

V2=U1*< 1 . r  A  a  (  2  > ) 

SnsF>d2  =  5n«5Pdl:!t  <  1  .  +Ar>  ( 3 )  ) 
T  2  =  T 1  It  <  1 .  +  A  a  (  4  >  ) 
Gasdn2=Gasdn  i  *  >  i  .  tAs  ( 5 ,'  ) 
P2=P1*< 1 . +Aa  <  o  > ) 


Re- in 1 1 i a  1 i ze  with  the  new  valu 

Msardl=Msard2 

U1=V2 

Sr.'»>*dt=;3  n  s  p  d  2 
T 1--T2 

G3sdnl=Gssdn2 
Mi:;l=U2 
P 1  =P2 

Rhowl =Rhow2 
Rhod3=G3sdnl-Rhowl 
El=RmwitRhowlitT  1 
Rn.=Rma/(  1 .-< 1  .H=>s)*El/Pl  > 
Cv=< 1 .+1 ,02*U2)*Cva 
Cp= ( 1 . +0 . 90*U2 ) *Cpa 
Gaml=Gam2 
Aresl=Area2 
Mdatl=Mdot2 
DO  93  1=1,20 
DAPLOT ( I ) =RA ( I ) *2. 0E6 
CONTINUE 
GO  TO  30 
i  CONTINUE 

IF  (NFL0T.LE.1)  GO  TO  101 

T1MAX=-1.E6 

V1MAX=-1.E6 

RHOMAX=-l . EA 

T1MIN=1 *E6 

U1MIN=1.E6 

RH0MIN=1 ♦ E6 

DO  350  1  =  1 » KPLOT 

T1MAX=AMAX1 (T1MAX,T1YYY<  I ) ) 

T1MIN=AMIN1 (T1MIN,T1YYY< I ) ) 

V1MAX=AMAX1 ( VI MAX  ,V1YYY(I) ) 

RH0MAX=AMAX1 ( RHOMAX ,  RHYYY < I ) ) 

U1MIN=AMIN1 <V1MIN,U1YYY< I ) ) 

RH0MIN=AMIN1 < RHGMIN , RHYYY < I ) ) 

'  CONTINUE 

INT1MX=INT ( T1MAX+0 .  5 ) 
IDEL=M0D( INTI MX, 5) 

T1MAX=INT 1MX+5-IDEL 
I NT 1HN=INT (T1MIN-0.5) 

IDEL=M0D ( INT1MN » 5  > 

T1MIN=T 1M IN-5+1 DEL 
INU1MX=INT <  V1MAX  f 0 . 5  > 
IDEL=MOD< I NO 1 MX  >5 ) 

U 1  MAX  =  I  NO  1  MX +5- 1  DEI. 

I NV1MN=INT (V1MIN-0.5) 
IDEL=MOD( INOlMNf 5> 
V1MIN=INVIMN-SHDEL 


IDEL=MOD\ INRHn  w  5 ‘ 

RHQMAX= INRHMX  S-5  [BEL 
INRHMN= INT <  RHDMIN-0 . 5  > 

IDEL=MOD< INRHMN.5) 

R  H  0  M I N  =  I N  R  H  M  N  -  5  +  I D  E L 
INPLOX=INT  <  PLTXXX  ( KP L 0 r  '■  f  o .  5 > 

IDEL-MOCK  INPLOX , 10) 

PLOXMX= INPLOX+ 1 0— I DEL 
XXSTP=PL0XMX/10. 

VI STP= ( '  ’  1MAX- V l MI N  > ,  1 V . 

T 1  STP=  <  T1MAX-  T 1 M 1 N  )  /  l. 0  - 
KOSTP=  (  RHOMAX-RHOM  i*,,  t-*. 

CALL  BGNPL ( KPLQT+ l ) 

CALL  PAGE  < 12. ,9. ) 

CALL  NOBRDR 

CALL  AREA2D< 10. >7.5) 

CALL  XNAME  ( ' DISTANCE  FROH  ENGINE , ( Meters )  * ' , 100 ) 
CALL  YNAME  ( ' TEMPERATURE , < Decrees  Cent israde ) t '» 100 ) 
CALL  GRAF  ( 0 . 0 » XXSTP , PLOXMX , T1MIN , T 1STP , T 1 MAX ) 

CALL  CURVE  < PLTXXX , T 1 YYY , KPLOT , 1 ) 

CALL  ENDPL  (KPLOT+1) 

CALL  BGNPL (KPLOT+2) 

CALL  PAGE  ( 12. f9. ) 

CALL  NOBRDR 

CALL  AREA2D  <10., 7.5) 

CALL  XNAME  ('DISTANCE  FROM  ENGINE *< Me ters )*', 100 ) 
CALL  YNAME  < 'EXHAUST  GAS  SFEEDr (Meters/Sec. )$' , 100) 
CALL  GRAF  ( 0 . 0 ,  XXSTP  ,  PLOXMX  ,  VI MIN ,  V.1STP  ,  V1MAX  ) 

CALL  CURVE  (PLTXXX , VI YYY , KPLOT . 1 ) 

CALL  ENDPL  (KPLOT+2) 

CALL  BGNPL  (KPLOT+3) 

CALL  PAGE  <12. ,9.) 

CALL  NOBRDR 

CALL  AREA2D  ( 10. ,7.5) 

CALL  XNAME  ('DISTANCE  FROM  ENGINE ,( Meters >*', 100 ) 


CALL  YNAME  ('GAS  DENSITY, (KS , /Cubic  Meter )*', 100 > 
CALL  GRAF  <0.0, XXSTP, PLOXMX, RHOMIN,ROSTP,RHOMAX> 
CALL  CURVE  (PLTXXX, RHYYY, KPLOT, 1 ) 

CALL  ENDPL (KPLOT+3) 

101  CONTINUE 

CALL  DONEPL 
202  FORMAT  (/) 

STOP 

END 

FUNCTION  Esat(T) 

C  DATA  Ts  /373.16/ 

C  DATA  D5  /3. 0057149/ 

C  D5=AL0G(EUS)  ,  where  EU3 =1013. 246  millibars. 

C  IF  (T.GT.Ts)  GO  TO  1 

C  A=T  s/T 

0  B=( 1 .-1 ./A>*11 .344 

C  C  =  - 3. 49149* < A -1 .  ) 

C  D1  =  -7.9O2?0*(A-1  .  ) 

C  D2  =  5 . 02808*AL0G 1 0  <  A ) 

C  D3--1 .3816E-7#10.#;l<fB-l  .  ) 

C  D4=8. 1328E-8*10.**<C~1  . 

C  Tt.=Dl  +D2+D3  +  D4  +  D5 
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C  E=10.*iTt 

i:  Ro  ro  2 

a  CONTINUE 

E=E'<P  ( 19  .  3450=2-4522  .  0246  •■'  r  ) 

2  CONTINUE 

Esat=100 . *L 
END 

FUNCTION  FNAEEA ( Z . mO ,  Rat  i  o ) 

ZZZZ-Z 
RA  TTT =Ra t 1 o 
FNAREA=AO 
END 

FUNCT  ION  ORACr  •  T-  •  '.'AIR  •  UDMJP  •  KHO  > 

UR-VAIR-VDKCP 

LiRAGF  =  3 . 129E-4*R*VR*  (  l  .  +2 . 249E4*RH0*R.*L'F  ) 

0  The  drag  force  is  calculated  fru.ii  Oarer,  flow  jpf  roxi .nation  • 

C  (draa  force)  =  6#PI*£ta#R*U#  <  1  K  3/ 1  6  )Re  )  »  where  Eta  is  the  v  t  •,  - 

C  cosity.R  is  the  radius  of  the  drop*  U  is  its  velocity  relative- 

C  to  the  air  flow?  Re  is  the  Reynold's  number »  defined  by 

C  Re=2#Rho*V/Eta .  where  Rho  is  the  density  of  the  air. 

END 

c  PROGRAM  MIESCAT 

c  DIMENSION  JAY ( 300 ) .UHY(300) , SIGH ( 300 >. CJAYPM ( 300 > .SIGHP(1(300) > 

C 

C  The  nomenclature  used  in  the  documentation  of  this  program 

C  follows  that  used  in  ‘HANDBOOK  OF  MATHEMATICAL  FUNCTIONS* .edited 

C  by  ABRAMOWITZ  and  STEGUN.  published  by  NATIONAL  BUREAU  OF  STAND- 

C  ARDS.  December. 1965. 

C 

C  JAY(N)  are  the  Spherical  Bessel  Functions  of  the  first  kind. 

C  UHY(N)  are  the  Spherical  Bessel  Functions  of  the  second  kind. 

C  SIGH<N)  are  the  Riccati-Bessel  Functions  (with  real  argument). 

C  CJAYPM(N)  are  the  first  derivatives  of  the  Spherical  Bessel 

C  Functions  of  the  first  kind. 

C  SIGHFM(N)  are  the  first  derivatives  of  the  Riccati-Eressel 

C  Functions. 

C 

C  1CPK300)  »CTAU(300)  »CDPI  (300) 

C 

C  CPI(N)  are  the  Nth  order  Legendre  Functions  divided  by  the 

C  Sine  of  the  scattering  angle. THETA  » 

C  CTAU(N)  are  the  first  derivatives  of  the  Legendre  Functions 

C  with  respect  to  the  scattering  angle  THETA. 

C  COPI(N)  are  the  first  derivatives  of  the  Legendre  Functions 

C  with  respect  to  the  Cosine  of  the  scattering  andle  THETA. 

C 

c  COMPLEX  CHI (300) . ETA ( 300 ) . CH1PRM ( 300 > . ETAPRM ( 300 ) .CSIGHR(300> . 

C 

C  CH1(N)  are  the  Spherical  Bessel  Functions  of  the  third  kind. 

C  ETA(N)  a.e  the  complex  Riccati-E<essel  Functions  with  complex 

C  arguments. 

C  CHIPRM(N)  3re  the  first  order  derivatives  of  the  Spherical  Hesse 

C  Functions  of  the  third  kind, 

C  ETAPRM(N)  are  the  first  order  derivatives  of  the  R iccati -Bessel 

C  Functions. 

C  C3IGHR(N)  are  the  ratios  of  SIGHPM(N)  to  SIGH(N)  For  comele  . 

C  arguments. 

C 
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c  Kj!  1  H * CS2TH  •  CCEXT  * CEM  •  CXn 

c 

CS1TH  is  the  S  sub  l  a  1. Client  of  th«i  amp".  itude  sc  sttenra  wh'i 
(  See  ‘ATMOSPHERIC  R  A  DIPT  ION*  •  v  r ,  1  •  l  •  b  s  R  .  M  .  GOODY  to-  the  defjr.i 
tion  and  nomenclature  of  the  scattering  eoust i on* .  ) 

CS2TH  is  the  S  sub  2  element  of  the  scattering  matrix* 

CCEXT  is  the  extinction  cross-section* 

CEM  is  tiie  complex  index  of  refraction  of  the  scstterer* 

L  CXM  is  tide  complex  wave  number  (  =-CEM*IGA  > 

(1 

c  REAL  JAY  * JAYO  *  KOA  *  KO  *  LAMBDA  *  MMM1 *  dMM2 

L 

KOA  is  the  wavenumber  ■.  *2-!tp  I /LAMBDA  >  > 

C  NO  is  2*P I /LAMBDA* 

C  LAMBDA  is  the  wavelength  < in  meters)  of  the  incident  radiation* 

C  MMM1  is  the  real  part  of  the  index  of  refraction* 

C  MMM2  is  the  imaginary  part  of  the  index  of  refraction, 

c  DATA  CPI i 1 ) t CDPI ( 1 ) rCDPI (2>  /-I . 0*0 . 0* -3 . 0/ 

c  DATA  PI *EPS0  rCSPLGT  /3. 1415927*  6 . 84 194 13E- 1 2  *  3.00E8  / 

c  DATA  RADIUS* LAMBDA  /  1.00E-6*  4123.0E-10  / 

c  DATA  MMM1 * MMM2  /  1.50  ,  -0.01  / 

c  DATA  THETA  /  0.000000  / 

c  DATA  RNDERR  /  1.00E-8  / 

c  K0=2.*P I /LAMBDA 

c  K0A=K0*RADIUS 

c  CEM=CMPL  X ( MMM1 *  MMM2 ) 

c  CXM=K0A*CEM 

c  TEST1=CABS(CXM) 

c  IF  ( (K0A.LT.0.5) .AND. (TEST1 .LT.0.5) )  GO  TO  50 

c  IF  (K0A.LE.1.0)  GO  TO  60 

c  RNDA=-AL0G( 10.** (-RNDERR ) ) 

c  RNDA323=(3.*RNDA>*(2./3. > 

c  B=RND323/10. 

c  NMAX=RND323 

c  IF  (KOA.GT.B)  NMAX= ( RND323+ < RND323**2+ ( 2 . 0*K0A ) **2 ) ** ( 0 . 5 ' ) /2 . 

c  UNEVNN-NMAX+O • 5 

c  TANHAL- ( 1 . - ( KOA/UNEVNN ) **2 ) ** ( 0 . 5 ) 

c  ALPHA=0 . 5*AL0G ( ( 1 . +TANHAL ) / ( 1 . +TANHAL ) ) 

c  JAY  ( NMAX )  =EXP  ( UNEUNN*  (  TANHAL-ALPHA ) ) / ( 2 . 0* (  K0A*TANHAL*UNEL'NN  ) 

c  1** ( 0 . 5 ) ) 
c  KMAX=NMAX-1 

c  IF  ( NMAX, GT. 299)  GO  TO  100 

c  JAY ( NMAX + 1 ) =0 ♦ 0 

c  SUMJ=(2*NMAX+1 >* JAY < NMAX )* JAY ( NMAX > 

c  DO  1  1=1* KM AX 

c  N2N1=2* ( NMAX-I+1 ) +1 

c  JAY ( NMAX-I >=N2N 1* JAY (NMAX-I+1 )/K0A-JAY( NMAX- I +2) 

c  SUM J=SUM J+ JAY ( NMAX-I ) *JAY i NMAX- 1 ) *N2N1 

c  1  CONTINUE 

C  JAY0=2 , 0* JAY ( 1 )/KQA-JAY (2) 

c  SUMJ=SUMJ+ JAYO* JAYO 

c  SUM J=SUM J** (0.5) 

c  JAY ( 1 ) = JAY ( 1 ) /SUM J 

c  WHY0= -COS ( KOA ) /KGA 

c  WHY  ( 1 )  =UHY O/KOA-S IN  (  KOA  '  i'  CA 

c  CH 1(1) =CMPLX ( JAY ( 1 > *  WHY . 1' ' 

c  WHY ( 2 ) =3 . *WH ( 1 ) /KOA  -WHYO 

c  C0STH=C0S( THETA) 
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JAYO  AND  JAY1 


c  SINTH=  <  i  .  -CUSTH-tCOSTH  :  **i  J.5> 

SIN  TH  =  S  I  i  J  (  THETA 
CPI  (2) -—3.  StCOS  T !  •: 
c  SINTH2=S INTHSSI.  <  IH 

c  CTAU< I >-~CQSTH 

c  DO  2  N=2»NMAX 

c  WHY  ( N+ 1 )  =  (  2*NT  1  )  YUHY  (  N  >  /KOA-WH  <  ( r>-  1) 

<:  JAY  (  N  )  -  JAY  <  N  )  /SUM  J 

c  CHI C  N ) -CMPLX ( JAY  v N ) , WHY <  N > ) 

c  CPI  (N+l )»(  (  2*N >  1  )  /N  )  ttCOS  TH-HCPI  (  N )  - (  C  N  +  l  )  Ji )  *CP I  \  ,\  -  !.  ) 

c  CTAU  ( N  )  =COSTH*CP  I  <  N  > -SIN  TM2  fcCDPI  ( N  > 

■'  CDPI(N+1)*(  (2KN-*-!)  '<N-t> )*CUSTH*CDPI<M  •  (  ( N+2 ) / < :4~  1  >  » 

C  J  CUNTLNCE 

c  CH 1 l NMA.<+ L ) -CMPLX (  >AY ( NMAX  f 1 ; • WHY ( NMAX+ i '  ) 

C  CHECK  ACCURACY  OF  JAYS  HERE  BY  JAYO  AND  JAY1 
c  CEM=CMPLX  < rtMM l »  MM.12  ’ 

c  CXM-KOA#CEM 

c  DO  3  N=1 » NhAX 

<?  SIGH<N)=KOA*JAY. N) 

c  ETA  (  N  )  =K0A*CH1  (  N ) 

c  C JAYPM (  N  )  =  (  N*  JA  f  <  N- 1 >  -(MH  >*JAY(N+1  )  )/(2*NH> 

c  CH1PRM(N)=(N*CH1 (N-l )-<N+i ) *CH1 (N+l ?  )/ « 2KN+1 > 

c  S1GHPM ( N ) = JAY ( N ) +KOA*C JA YPh ( N  > 

c  ETAPRM  ( N )  =CH  1  (  N )  +K0A*CH1PRM  (  N  > 

c  CSIGHR <.  N ) =FMCALC ( CXM  » N *  RNDERR ) +N/CXM 

C  3  CONTINUE 
c  CS1TH=<0.0»0.0) 

c  CS2TH= ( 0  4  0 » 0 . 0  > 

c  CCEXT  =  <  0 . 0  r  0 . 0 ) 

<:  DO  4  N~1  >  NhAX 

c  NFAC2 1=2*N+ 1 

c  NFAC12=NFAC21/(N*(N+1) > 

c  CS1TH=CS1TH4  NFAC12*(CA,1IE(H)1!CPI  ( N ) +CBMIE ( N > *CTAU ( N )  ) 

c  CS2TH=CS2TH+NFAC1 2*  (  CAMIE  ( N  )  *CTA»J (  N  >  +CBM IE  (  N )  *CP I  ( N  )  ) 

c  CCEXT=CCEXT+(2*N+1)*(CAMIE(N>+CBMIE(N) > 

c  QSCA=QSCA+NFAC21*  ( CAMIE  (  N  )  *C0N JG  < CAMIE <  N  )  + 

c  1CBMIE(N>*C0NJG(CBMIE<N) ) 
c  4  CONTINUE 

c  CS1TH=CS1TH#ELEC0*<0,«0*-1 .0) 

c  CS2TH=CS2TH*EL£C0*<0.0r-l .0) 

c  CEXT=(2.*PI/(K0*K0) )*REAL( CCEXT) 

c  0EXT=(2./(K0A*K0A> ) *REAL < CCEXT ) 

e  QSCA=QSCA* ( 2 . / ( KOA*KOA ) ) 

c  QABS=GEXT-GSCA 

cccPUT  IN  PRINT  STATEMENTS  HERE 
c  GO  TO  200 

c  50  CONTINUE 

c  CMPLF1-  ( CEM:*CEM-1 .  ) / ( CEM*CEM+2 .  > 

c  CS1TH=KQA**3*CMFLF1 

c  CS2TH=CS1TH*C0S( THETA) 

e  REEL=CABS ( CMPLF 1 ) 

c  QSCA-(3./3.  )XK0A**4*REEUtREEL 

c  GAeS=»4  .  KAIMAG  (  CMPLF  1 ) 

CCC  PRINT  out  RAYLEIGH  scatterma  here, 
c  GO  TO  200 

■'  60  CONTINUE 

c  A=MMM1 

C  E«=MMM2**2 


Zt=<A+B>**2+4.*< A-B>+4 . 

Z2 -4 .  *  ( A+B )  **2v  1 2  .  %  '  A-B  >  r  ?  . 

QSCA=3 .*< ( ( A+B  > **2 f A-B-Z . ) *42 +36 . *A*B > *K0AS*4 
l*<  1 .  +  < 1 .2/2  ■  *  ;x  <  A+B  -4 .  ) *K0A**2-8 . *MMMl*MMM2*KnA  5*3/21 )  /■ 
QEXT=  ( 24  ♦  *MMH 1*MMM2*K0 A/Z i )  +  <  *  4  •  /  '5 .  >  +  <20 ./ ( 3 .  *Z2  >  '  M3: 
t  ( 7 . 4  ( A+B ) **2+4  *  4  <  A— B— 5 ♦ > ) ) *MMM1*MMM2*K0A*:43 
2+<8./(3.*Zl*42> )4( (  <A+B>**2+A-B-2.  >**2+36.*A#B>*itOA**4 


,*Z’.* 


C  ^-SS^QEXT-GSCA 

CCCCCCC  P-TNT  OUT  HIGHER  ORDER  EXPANSION  RESULTS  HERE 
c  GO  TO  200 


c  100  CONTINUE 

ccccPRINT  'NUMBER  OF  TERMS  REQI'~  "~D  TO  CALCULATE  THE  MTE  COEFG.TCO 
c  200  CONTINUE 

c  COMPLEX  FUNCTION  FrtCALC ( X » N f RNDERR ; 

c  COMPLEX  ClXf X»CNUH(300) > COEN < 300 >  » CAM  1 300  > * FhCALC  * CNUF 

c  r  '>  CDEN(l)  /  (1.0»0.0>/ 

c  CIX=C  M/X 

c  DO  1  M=1 »300 

c  CAM(M>=-<2*(N+M>  *CIX 

c  1  CONTINUE 
c  CNUM ( 1 ) =C  AM  < 1 ) 

c  CNUM ( 2 ) “CAM  <  2  >  + 1  • /CAM  < 1 ) 


c  CDEN(2)=CAM(2: 

c  FMCALC -CNUM ( 1 ) 4 CNUM C  2 ) /CDEN ( 2 ) 

c  DO  2  M=3»300 

c  CNUM  <  M  )  =CAH  (  M  )  -f  t  .  /CNUM  <  M- 1 ) 

c  C DEN « M ) =CAM  <  M  >  +  1 ♦ /CDEN <  M- 1  * 

c  CNUF -CNUM ( M ) /CDEN  <  M ) 

c  FMCALC=FMCALC  SC'fUF 

c  F«CABS<CNUF> 

c  AB3F=AF:S<F-1  .0 

IF  «  A3SF  *LT  .P>"  f-RR)  GO  TO  3 
c  2  CONTINUE 
c  3  COMTINL"" 
c  RETURN 

c  END 


A-?1 


